Let us consider the dataframe wages, which containes information about 3294 USA working individuals. The data are taken from the National Longitudinal Survey and are related to 1987. The variable as are listed below and the output of the str command is given
exper
experience in years
male
1 male, 0 female
school
years of schooling
wage
wage (in 1980$) per hour
region
Center, North, South
## 'data.frame': 3296 obs. of 5 variables:
## $ exper : int 9 12 11 9 8 9 8 10 12 7 ...
## $ male : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...
## $ school: int 13 12 11 14 14 14 12 12 10 12 ...
## $ wage : num 6.32 5.48 3.64 4.59 2.42 ...
## $ region: Factor w/ 3 levels "Center","North",..: 1 1 3 1 1 1 1 1 1 3 ...The aim of the study is to analyze the potential relationship between the response variable wage and the explanatory variables considered in the dataframe. Describe how to perform a preliminary data analysis on this dataframe, using suitable R commands. Moreover, consider the following plots and discuss the possibility of measuring the variable wage in the logarithmic scale
Distribuzione delle variabili
source("./../functions.R", local = knitr::knit_global())
summary(wages)## exper male school wage
## Min. : 1.000 female:1569 Min. : 3.00 Min. : 0.07656
## 1st Qu.: 7.000 male :1727 1st Qu.:11.00 1st Qu.: 3.62673
## Median : 8.000 Median :12.00 Median : 5.20578
## Mean : 8.042 Mean :11.63 Mean : 5.81639
## 3rd Qu.: 9.000 3rd Qu.:12.00 3rd Qu.: 7.30723
## Max. :18.000 Max. :16.00 Max. :112.79193
## region
## Center:1164
## North :1105
## South :1027
##
##
## numeric.graph(wages, c(4,1,3), qqplot = 4, no.density = c(1,3), draw.normal.curve = c(1,3))
## [1] "Skewness indexes"
## wage exper school
## 9.5491015 0.1555681 -0.3934935
## [1] "Kurtosis indexes - 3"
## wage exper school
## 212.9197746 0.6528604 1.2394911categorical.graph(wages, c(2,5))
La variabile risposta presenta una elevatissima simmetria e curtosi, con i valori praticamente tutti concentrati a sinistra del grafico con la presenza di due outlier elevatissimi. Il qqplot conferma la presenza di una coda destra molto pesante, mentre i quantili centrali seguono abbastanza bene la distribuzione normale. Si consiglia perciò o di rimuovere gli outlier più lontani oppure di introdurre una trasformata.
exper invece sembra seguire piuttosto bene la curva della distribuzione normale (buona simmetria, curtosi leggermente ipernormale).
Considerazioni simili valgono anche per school (anche qui buona simmetria), la distribuzione di presenta con una frequenza del valore centrale (11) molto elevata e una coda sinistra leggermente appesantita.
Vista la natura di queste due variabili non ha molto senso considerare gli outlier.
wages$logwage <- log(wages$wage)
numeric.graph(wages, 6, qqplot = 6)
## [1] "Skewness indexes"
## [1] -1.11844
## [1] "Kurtosis indexes - 3"
## [1] 4.28369La trasformata logaritmica migliora notevolmente sia simmetria che curtosi della variabile, tuttavia appesantisce molto la coda sinistra (infatti la campana rimane fortemente ipernormale), per cui eventuali assunzioni di normalità sulla variabile vanno fatte con cautela.
Relazioni tra variabili
boxplots.graph(6, c(1:3,5), wages)
Dai grafici non sembrano esserci delle particolari correlazioni tra log(wage) e le altre variabili.
boxplots.graph(1, 2, wages)
mean(wages$exper[wages$male == 'male'])## [1] 8.323104mean(wages$exper[wages$male == 'female'])## [1] 7.732314Dal grafico e dal confronto delle medie sembra esserci una (debole) correlazione tra exper e male.
Introduzione di trasformate rimuovendo gli outlier
#Rimuovo outlier
newWage <- wages[wages$wage < 30,]
newWage$logwage <- log(newWage$wage)
newWage$wageT <- transformResponseWithBoxCox(newWage$wage~1, newWage, 4)
## [1] "Ideal lambda for wage : 0.4"numeric.graph(newWage, c(6,7), qqplot = c(6,7))
## [1] "Skewness indexes"
## logwage wageT
## -1.23982911 0.04971715
## [1] "Kurtosis indexes - 3"
## logwage wageT
## 4.1493448 0.8880249La trasformata di Box Cox con \(\labda=0.4\) sembra preferibile rispetto a quella logaritmica.
In order to describe the effect of the factor male on the response log(wage) we may analyze this plot, where the probability distribution of log(wage) is represented by considering the kernel density estimates conditioned on the two levels (1 male, 0 female) of the variable male
Il grafico conferma l’assenza di correlazione tra log(wage) e male, infatti le due curve di densità sono praticamente sovrapponibili.
With the commands mod.0<-lm(log(wage) ∼ male,data=wages) and mod.1<-lm(log(wage) ∼ exper*male, data=wages), two regression models are defined for describing the potential effect of male and exper on the response log(wage). Comment the model fitting outcomes given by the function summary (Hint: consider the fact that the average years of experience in the sample is lower for women than for men).
mod.0 <- lm(log(wage) ~ male,data=wages)
summary(mod.0)##
## Call:
## lm(formula = log(wage) ~ male, data = wages)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0445 -0.3073 0.0544 0.3839 3.0325
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.47475 0.01559 94.59 <2e-16 ***
## malemale 0.21826 0.02154 10.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6176 on 3294 degrees of freedom
## Multiple R-squared: 0.03023, Adjusted R-squared: 0.02994
## F-statistic: 102.7 on 1 and 3294 DF, p-value: < 2.2e-16mod.1 <- lm(log(wage) ~ exper*male, data=wages)
summary(mod.1)##
## Call:
## lm(formula = log(wage) ~ exper * male, data = wages)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0906 -0.3050 0.0560 0.3792 3.0468
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.193551 0.057470 20.768 < 2e-16 ***
## exper 0.036367 0.007156 5.082 3.94e-07 ***
## malemale 0.463707 0.079062 5.865 4.93e-09 ***
## exper:malemale -0.032071 0.009518 -3.369 0.000762 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6153 on 3292 degrees of freedom
## Multiple R-squared: 0.03792, Adjusted R-squared: 0.03704
## F-statistic: 43.25 on 3 and 3292 DF, p-value: < 2.2e-16Il confronto tra i due modelli mostra come l’interazione tra exper e male (rappresentato dalla variabile dummy malemale) sia rilevante nel modello, di conseguenza lo sono anche i loro main effects.
Finally, a complete regression model is fitted mod.2 <- lm(log(wage) ∼., data=wages) and the following output is obtained by the R function summary.
summary(lm(log(wage) ~. -logwage, data=wages))##
## Call:
## lm(formula = log(wage) ~ . - logwage, data = wages)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0008 -0.2821 0.0468 0.3673 3.2337
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.279709 0.090321 -3.097 0.00197 **
## exper 0.034618 0.004549 7.610 3.55e-14 ***
## malemale 0.246474 0.020607 11.961 < 2e-16 ***
## school 0.122909 0.006278 19.578 < 2e-16 ***
## regionNorth 0.051107 0.024505 2.086 0.03709 *
## regionSouth 0.047168 0.024969 1.889 0.05898 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5831 on 3290 degrees of freedom
## Multiple R-squared: 0.1364, Adjusted R-squared: 0.1351
## F-statistic: 103.9 on 5 and 3290 DF, p-value: < 2.2e-16Describe how to interpret these results, and then suggest how to proceed with further analyses.
Il modello qui presentato mostra come tutti le variabili siano regressori rilevanti per logwage.
plot(lm(log(wage) ~. -logwage, data = wages))
Il modello presenta alcune criticità dal punto di vista della normalità dei residui. Le ipotesi di omoschedasticità e linearità sono rispettate.
Modello alternativo
Cosa otteniamo utilizzando sempre la trasformata logaritmica senza outliers
fitLog <- lm(logwage~exper+male+school+region, data = newWage)
summary(fitLog)##
## Call:
## lm(formula = logwage ~ exper + male + school + region, data = newWage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9969 -0.2792 0.0490 0.3710 1.8645
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.299749 0.089186 -3.361 0.000786 ***
## exper 0.036006 0.004493 8.013 1.54e-15 ***
## malemale 0.243299 0.020352 11.954 < 2e-16 ***
## school 0.123402 0.006196 19.915 < 2e-16 ***
## regionNorth 0.056040 0.024193 2.316 0.020597 *
## regionSouth 0.045451 0.024669 1.842 0.065500 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5754 on 3285 degrees of freedom
## Multiple R-squared: 0.1405, Adjusted R-squared: 0.1392
## F-statistic: 107.4 on 5 and 3285 DF, p-value: < 2.2e-16plot(fitLog)
drop1(lm(wageT~exper+male+school+region, data = newWage), test = "F")## Single term deletions
##
## Model:
## wageT ~ exper + male + school + region
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 3300.7 21.66
## exper 1 55.62 3356.3 74.66 55.3587 1.275e-13 ***
## male 1 170.70 3471.4 185.60 169.8906 < 2.2e-16 ***
## school 1 444.86 3745.5 435.76 442.7430 < 2.2e-16 ***
## region 2 3.29 3304.0 20.94 1.6363 0.1949
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dal risultato di drop1 risulta che per la nuova trasformata il regressore region non è significativo.
fitBoxCox <- lm(wageT~exper+male+school, data=newWage)
summary(fitBoxCox)##
## Call:
## lm(formula = wageT ~ exper + male + school, data = newWage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7995 -0.6243 -0.0142 0.6259 4.4220
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.997868 0.153070 -6.519 8.16e-11 ***
## exper 0.057779 0.007824 7.384 1.93e-13 ***
## malemale 0.462388 0.035457 13.041 < 2e-16 ***
## school 0.226935 0.010796 21.021 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.003 on 3287 degrees of freedom
## Multiple R-squared: 0.1523, Adjusted R-squared: 0.1516
## F-statistic: 196.9 on 3 and 3287 DF, p-value: < 2.2e-16plot(fitBoxCox)
La trasformata calcolata con il metodo di Box-Cox permette di rientrare (seppur parzialmente) nell’ipotesi di normalità dei residui.
AIC(fitLog, fitBoxCox)## df AIC
## fitLog 7 5710.026
## fitBoxCox 5 9362.390BIC(fitLog, fitBoxCox)## df BIC
## fitLog 7 5752.719
## fitBoxCox 5 9392.884Gli indici \(R^2\) sono più alti nella trasformata Box Cox rispetto al logaritmo, tuttavia le statistiche AIC e BIC tendono a preferire quest’ultima.